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, Abstract. Digital deblurring of images is an important problem that arises in multifrequency observations of 

the Cosmic Microwave Background (CMB) where, because of the width of the point spread functions (PSF), 
maps at different frequencies suffer a different loss of spatial resolution. Deblurring is useful for various reasons: 
. first, it helps to restore high frequency components lost through the smoothing effect of the instrument's PSF; 

CL(' second, emissions at various frequencies observed with different resolutions can be better studied on a comparable 

JL I resolution; third, some map-based component separation algorithms require maps with similar level of degradation. 

■ Because of computational efficiency, deblurring is usually done in the frequency domain. But this approach has 

some limitations as it requires spatial invariance of the PSF, stationarity of the noise, and is not flexible in 
, the selection of more appropriate boundary conditions. Deblurring in real space is more flexible but usually not 

used because of its high computational cost. In this paper (the first in a series on the subject) we present new 
algorithms that allow the use of real space deblurring techniques even for very large images. In particular, we 
, consider the use of Tikhonov deblurring of noisy maps with applications to PLANCK. We provide details for 

5_( ' efficient implementations of the algorithms. Their performance is tested on Gaussian and non-Gaussian simulated 

CMB maps, and PSFs with both circular and elliptical symmetry. Matlab code is made available. 

Key words. Methods: data analysis - Methods: statistical - Cosmology: cosmic microwave background 



1. Introduction 

During the last decade, observations of the Cosmic 
Microwave Background (CMB) anisotropics have pro- 
gressed significantly. After the first evidence of CMB in- 
tensity fliictuations measured by the COBE satellite (see 
Smoot I I1999I and references therein) , several balloon- 
borne and ground based experiments have successfully de- 
tected CMB anisotropics at degre e and sub-degree aiigular 
scales fae Bernardis et al. I l2002t iHalverson et al. I l2002t 
[Lee et al. ..2001: Padin et al. Il200l[l . The MAP satellite ^ 



Send offprint requests to: R. Vio 
^ http://map.gsfc.nasa.gov 



currently in operation will soon provide full sky maps of 
CMB anisotropics at about 20 arcmin. resolution and a 
sensitivity of the order of 10 /iK, on a frequency range 
extending from 22 to 90 GHz. This extraordinary ex- 
perimental enterprise will be followed by the PLANCK 
satellite, scheduled for launch in 2007 ^ , that will provide 
full sky maps of total intensity and polarization of CMB 
anisotropy at a few arcmin. resolution, and a sensitivity 
of a few /LtK, on nine frequency channels ranging from 30 
to 857 GHz. 

These new sets of observations will pose new and chal- 
lenging questions in data analysis. Methods will have to be 

^ http:/ /astro. estec.esa.nl/SA-general/Projects/Planck 
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developed to process the large amount of incoming data, 
and to extract and separate cosmological information from 
foreground emissions from extra-Galactic sources as well 
as from our own Galaxy. 

An important question is the creation of sky maps 
based on small angular scale time ordered data from 
all-sky CMB observations of MAP and PLANCK; effi- 
cient map-making algorithms based o n a maximum like- 
lihood approach have been proposed llNatoli et al. iboOlt 
iBorrill et aLlbnOlt IStomnor et al. II2OO2I) . Reefarding the 
separation of emissions coming from different astrophys- 
ical processes, multifrequency observations can be ex- 
ploited. For example, available prior information about 
the signals can be used in a regularised inversion via 
Wiener filtering and maxi mum entropy m ethods, either on 
smal l sky patches (.Hobson et al. lll998l) or on the wh ole 
sky l)Bouchet fc Gispert Ill999l: IStolvarov et al. Il2002() . It 
has been shown that under certain independence assump- 
tion on the signals, the map-operating algorithms based on 
Independent Component An alysis (ICA) techniques can 
be applied on sky patches teaccigalupi e t al~l l200(ll) as 
well as on the whole sky (^Maino et al 200^ 

In this paper we study a difi^erent problem. We consider 
the effects of the beam of the instrument used to gather 
the observations. We present efficient numerical methods 
to estimate the emission pattern lost through the degrada- 
tion of the instrument's point spread function (PSF) and 
noise contamination. This "deblurring" process may prove 
very useful in CMB data analyses: first, it helps recover 
high frequencies smoothed out by the instrument's PSF. 
Second, a better understanding of sky emissions, from 
foregrounds in particular, is achieved if multifrequency sky 
maps are compared on a common resolution. Third, some 
map-based component separation algorithms, such as ICA 
I Baccigaluni et al. I lionct iMaino et 111120021) require in- 
put maps with similar level of degradation. Furthermore, 
although the aim of satellite missions such as Planck and 
MAP is to obtain full sky maps of the CMB, the strength 
of the CMB over other backgrounds or contaminating 
sources will vary over the sky. Therefore, even if some 
characteristics of CMB are estimated on full sky maps, 
it will be convenient to check these results on smaller 
sky patches where CMB dominates the other components 
(e.g., at high Galactic latitude and at high observing fre- 
quency) and data are free from instrumental and/or ob- 
servational problems. 

Most deblurrin g of CMB data has been done in fre- 
quency space, fsee iHobson et al."1ll998t Fstolvarov et al. I 
[2002.1 .' This approach is computationally efficient but has 
some important limitations; it requires the stationarity of 
the contaminating noise and the spatial invariance of the 
PSF. Furthermore, it implicitly assumes periodic bound- 
ary conditions which, as we show below, is not necessarily 
the best choice. Tikhonov regularization provides a more 
flexible real domain deblurring technique that can be used 
when any of these assumptions are not met. 

Until recently, Tikhonov deblurring was avoided be- 
cause of its computational cost. For example, for a.n NxN 



pixel image in the frequency domain one works with NxN 
matrices, while in the spatial domain the matrices are of 
dimension N'^xN'^. However, new efficient algorithms that 
overcome this problem have made Tikhonov deblurring a 
competitive alternative to frequency domain techniques. 
In this first paper, we present some of these algorithms and 
show their good performance not only in regards to their 
computational cost but also in regards to their numeri- 
cal stability, which is an important characteristics in any 
ill-posed inversion problem. Furthermore, we will present 
some preliminary results concerning the deblurring for the 
case of spatially varying PSFs. A detailed treatment of this 
problem will be provided in a future paper. 

The rest of the paper is organized as follows. In Sects. El 
and 01 we provide a formal definition of the map-based 
deblurring process, focusing on boundary conditions and 
numerical issues. In Sect. 01 we discuss a regularization 
approach that is efficient and fiexible for deblurring noisy 
maps; results of numerical experiments are provided in 
Sect.ini Our first appfication to a reahstic map is presented 
in Sect. El w here we consider a r egion in the sky (already 
addressed bv lBaccigalupi et al. ti2000,.'l in which the CMB 
emission largely dominates over foregrounds. We will as- 
sume observational conditions, frequencies, angular reso- 
lution and noise level corresponding to the Low Frequency 
Instrument (LFI) of the PLANCK satellite. Further sim- 
ulations are presented in Sects. [71 and where we consider 
non- Gaussian random fields and spatially varying PSFs. 
In Sect. Owe close with final comments and conclusions. 

2. Formalization of the problem 

When a two-dimensional object /(C,??) is observed 
through an optical (linear) system, it is seen as an im- 
age g{x,y) 



+ 00 -\-OG 



(1) 



where the function h(x, y, ^, 77), called a point-spread func- 
tion (PSF), represents the blurring action of the optical 
instrument. Model Q]) allows the PSF to vary with posi- 
tion in both [x, y) and (^, 77) variables [space- variant PSF). 
Often, it is possible to simplify this model by assuming 
that the PSF is independent of position [space-invariant 
PSF), and thus 



(-00 +00 



g{x,y) 



h{x-U-v)f{Ul) d^d7^. (2) 



Another useful simplification occurs when PSF is separa- 
ble, which means that it can be written in the form 

Hx,y,£„T]) = hi{x,£,) h2{y,'ii), (3) 

for the space-variant PSF, and 

h{x-£„y~T])^hi{x~£,) h2{y~T]), (4) 
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for the space-invariant PSF. fn this case models (Q), 
can be simplified to 



+00 

gix,y) = J hi{x,^) 

— 00 

and 

+ CXD 

g{x,y) = / hi{x - 



+ 00 



dC, (5) 



+ CXD 



h2{y ~v)f{^,v) dv 



dC, (6) 



respectively. We see that separability of the PSF implies 
independent blurring along the horizontal and vertical di- 
rections. 

Models (P), I©, and © are only theoretical, fn 
practical applications, we only have discrete noisy obser- 
vations of the image g{x, y), which we model as a discrete 
linear system 

g = Hf + z. (7) 

Here, g = vec(G) and / = vec(F) are one-dimensional, 
column arrays containing, respectively, the observed im- 
age G and the true images F in stacked order z is an 
array containing the noise contribution (assumed to be 
of additive type), and H is a, matrix that represents the 
discretized blurring operator. 

By image deblurring we mean solving for / given the 
linear system ((JJ. This process is not trivial for the fol- 
lowing reasons: 

1. Images are available only in a finite bounded region. 
However, points near the boundary are affected by 
data outside the field of view. 

2. If the observed image is an iV x iV matrix, then /, g, 
and z are N'^ x I vectors and H is an N'^ x N'^ matrix. 
This means that even for small images (e.g., N = 512), 
the matrix H can be quite large. 

3. The models ©, and © are particular ex- 
amples of a class of i ntegral equations whic h are well 
known to be ill-posed l|Wing fc Zahrt Ill99j) . This im- 
plies that the matrix H is severely ill-conditioned and 
standard techniques that do not take this into account 
are likely to fail. 

Points (1) and (2) are addressed in Sect.O and point 
(3) in Sect. El 

3. Numerical issues in image deblurring 

3.1. Boundary conditions 

To resolve the difficulty of points near the image boundary 
being affected by data outside the field of view, we have 
to impose some boundary conditions. In image processing, 
one of the following three boundary conditions is (either 
explicitly or implicitly) typically made: 



• Periodic boundary conditions imply that the image re- 
peats in all directions. That is, we assume the image 
X has been extracted from a larger image that looks 
hke: 

XXX 
XXX 
XXX 

• Zero boundary conditions imply that the scene outside 
the borders of the image X are all zero. That is, we 
assume the image X has been extracted from a larger 
image of the form: 


0X0 


• Reflexive boundary conditions model the scene outside 
the image boundaries as a mirror image of the scene 
inside the boundaries. That is, we assume the image 
X has been extracted from a larger image that looks 
hke: 

XfQ Xj. Xj.Q 
Xc X Xc 

XfQ Xrp XfQ 

where Xc is obtained by "flipping" the columns of X, 
Xr is obtained by "flipping" the rows of X, and Xrc 
is obtained by "flipping" the rows and columns of X . 

For a spatially invariant PSF, each of these choices im- 
poses a particular kind of structure on the matrix H . To 
describe these structures we need the following notation: 

• An N X N matrix with entries that are constant on 
each diagonal is called a Toeplitz matrix, and it can be 
written as: 



(8) 



\ /jtv-i hN-2 • • • ho J 

li aii N X N Toeplitz matrix has the additional prop- 
erty that each column (and row) is a circular shift of 
its previous column (row), then it is called a circulant 
matrix and it has the form: 

hi ho /lAT-i 
h2 hi ho 



( ho 


h-i 


h-2 ■ 


■■ hi^ 


N 


hi 


ho 


h-i ■ 


•• h2- 


-N 


h2 


hi 


ho ■ 


■■ hs- 


-N 



hi\ 

h2 
h3 



(9) 



\ /lAT-i hN-3 ■ ■ ■ ho / 

An N X N matrix with entries that are constant on 
each anti-diagonal is called a Hankel matrix. It can be 
written as: 

(ho hi h2 ■■■ hN-i \ 
hi h2 hs 
/i2 hs hi 



^ We recall that for a A'' x M matrix P, vec(P) 
(pf P2 ■ ■ ■ pIi)'^ with the i-th column of matrix P. 



hw 
hw+i 



(10) 



h2N-2 j 
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In two-dimensional applications, such as image deblurring, 
the matrix H is usually represented in block form. In this 
case, we need to specify the block structure of the matrix, 
as well as the structure of each block. For example, if a 
matrix has the form: 



H 



T — 



( Ho 


Hi 


H 2 


Hi 


n\ 


Hi 


Ha 


Hi 


H2 


N 


H2 


Hi 


Ho 


H3 


N 



N-l 



H 



N-2 



H 



(11) 



N-3 



Ho I 



where each {Hi} is an N y. N matrix, then we say Ht 
is a block Toeplitz matrix. If, in addition, each {Hi} is 
itself a Toeplitz matrix, then we say Ht is a block Toeplitz 
matrix with Toeplitz blocks. Note that we can combine any 
number of structures. In this paper, we need the following 
combinations: 

BTTB (block Toeplitz with Toeplitz blocks). 
BCCB (block circulant with circulant blocks). 
BHHB (block Hankel with Hankel blocks). 
BTHB (block Toeplitz with Hankel blocks). 
BHTB (block Hankel with Tocphtz blocks). 

With this notation, we can precisely describe the structure 
of H for each of the boundary conditions: 

• With periodic boundary conditions H is BCCB. 

• With zero boundary conditions H is BTTB. 

• With reflexive boundary conditions ff is a sum of 
BTTB, BTHB, BHTB and BHHB matrices. (Each of 
these matrices takes into account, to some extent, con- 
tributions from X, Xr, Xc, and X^cj respectively.) 

Very often the choice of boundary conditions is considered 
secondary to the development of an efficient deblurring al- 
gorithm. In fact, in many situations the type of boundary 
conditions is implicitly imposed by the algorithm itself. 
For example, Fourier methods implicitly assume periodic 
boundary conditions. However, an improper selection of 
boundary conditions may introduce edge effects such as 
discontinuities, which appear as a slow decay of Fourier 
coefficients and Gibbs oscillations at the jump, which has 
obvious consequences on nonGaussianity tests and reliable 
estimations of the regularization parameter (see below). 
Consequently, the choice of the boundary conditions must 
be considered an integral part of the deblurring operation 
and not merely a byproduct of some algorithm. 

3.2. Numerical methods for the solution of large linear 
systems 

Even for a very large matrix H, there are characteristics 
of this matrix that make finding a solution of the system 
Q) feasible. For example: 

• if is a highly structured matrix, then it is not nec- 
essary to store all its entries. For instance, for a BTTB 
matrix it is sufficient to store the the first column and 



first row, whereas for a BCCB matrix it is enough to 
keep the first row. 
• significant memory savings can also be achieved if most 
of the entries of H are zero, i.e. H is a sparse matrix. 

Both of these properties can be efficiently exploited to 
reduce the computational cost considerably. Further dis- 
cussion of these issues can be found in Appendix IbI 



4. Numerical regularization 

From Eq. it is clear the the instrument's PSF smooths 
out high frequency components of the signal, and thus 
high frequency information is lost. An important conse- 
quence of this is the lack of a unique solution of the linear 
system {TJ; any solution subjected to high frequency per- 
turbations will fit the data g equally well. This makes 
the problem of recovering the signal / ill-posed. The ill- 
posedness also affects the stability of the solutions; a small 
perturbation of the data may result in a completely differ- 
ent solution. Mathematically, the system is ill-posed when 
the singular values of H decay to zero too fast. The larger 
the smoothing effect of the PSF the faster the decay of 
the singular values. Small singular values lead to solutions 
which fit the data but have very large energy; a property 
that is not easily justified (for ot her illustrative examples 
of ill-posedness sec Hansen 1997! lTenorio I2OOII) . To find a 
stable meaningful solution we have to use a regularization 
method. 

Ill-posed problem can be regularized by imposing con- 
straints on the unknown signal. For example, a bound on 
the energy or a smoothness constraint. This class of con- 
straints can be easily implemented in the classical frame- 
work of Tikhonov regularization ijTikhonov Sz Arsenin I 
Il977|) . The idea is to find a solution f^^ that minimizes a 
weighted combination of data misfit (residual norm) and 
signal constraint: 



fx = argmin ( |j Hf 



-aWl 



(12) 



where A is a scalar quantity and L is, for example, the 
identity matrix (for energy bound) or a discrete derivative 
operator of some order (for smoothness constraint). The 
optimal smoothing functional for dcconvolution operators 
on an in finite domain w as determined by Aref'eva (1971 
(see also lCullum Ill97i^ who showed that it is completely 
determined by the decay rate of the Fourier coefficients of 
the unknown function. Unfortunately, in real applications 
this result cannot be used directly because, in addition to 
data being available only on a finite domain, it requires 
the function one is trying to recover. 

In two-dimensional image processing applications, L is 
often chosen to be a discretized Laplace operator, which 
is defined as the second order partial derivative 



V/(x,2/) 



9V 



(13) 
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A discrete approximation c an be i mplemented as a convo- 
lution of / with the kernel ijjain 1^989, p.353) 



(14) 



The implementation of this convolution and the precise 
form of the matrix L depend on the boundary conditions. 
In general, we have 




/Ti -I 

-I T2 



\T3 



■■ I 

/ T2 / 



(15) 



where I is the identity matrix, and the matrices Ti, T2 
and r3 depend of the boundary conditions: 

• for periodic boundary conditions T3 = —J, 

/ 4-1 -1\ 



Ti 



-1 



■. -1 
-1 4/ 



(16) 



• for zero boundary conditions T3 = 0, 

/ 4 -1 



= -1-2 = 



-1 



V 



•. -1 

-1 4 / 



(17) 



for reflexive boundary conditions T3 = 0, 



/ 



■. -1 

1 3 -1 

-1 



(18) 



2/ 



-1 

4 -1 
-1 ■•• 



-1 



V 



-1 

4 -1 
-1 3; 



(19) 



Other oper ators, such a s first order derivatives, can be 
used for L ijjahne Ill997j) . However, these operators have 
the disadvantage that they are typically applied in only 
one direction (e.g., x-derivative or y-derivative), and 



therefore are not isotropic. Nonlinear isotropic implemen- 
tations ar e possible bu t they are computationally more 
expensive dJahne Ill997l) . In contrast, the Laplacian is an 
isotropic linear operator that can be efficiently incorpo- 
rated into our algorithms (see Appendix El- 

In Eq. (|12|l . the regularization parameter X controls the 
weight given to the signal constraint relative to the data 
misfit. A good selection of A is critical. A large A favours 
a small solution norm at the cost of a large data misfit, 
while a small A leads to data overfit. 

We next consider practical implementations to both 
solve problem H12I) and estimate an "optimal" value of A. 

4.1. Numerical methods for Tikhonov map deblurring 

As previously mentioned, the image deblurring problem 
is severely ill-conditioned and regularization is needed in 
order to compute solutions that are not completely cor- 
rupted by noise. One approach is Tikhonov regularization, 
where the solution f^^ is given by Eq. (fT^ . Stable algo- 
rithms for computing this solution are usually developed 
by reformulating (|12|l as the damped least squares prob- 
lem 



mm 



H 

XL 



f 



(20) 



It is important to note that the Tikhonov method is just 
one approach that can be used for regularization but many 
other schemes may be used. For large scale image deblur- 
ring problems, iterative regularization methods, such as 
conjugate gradients, Landweber iteration, or expectation- 
maximization (sometimes referred to as Richardson- 
Lucy), are often recommended. Regularization is enforced 
through iteration truncation; that is, the iteration index 
acts as the regularization parameter. Although the specific 
details of various iterative methods may be different, they 
usually have the common property that the most compu- 
tationally expensive operation at each iteration is matrix 
vector multiplications with H. For large scale image de- 
blurring problems, these computations can be done very 
efficiently using fast Fourier transforms. 

A disadvantage with using iterative methods is that it 
is very difficult to determine when to stop the iteration; 
that is, how to choose a good regularization parameter. 
Many methods have been proposed, but they are usually 
most reliable when used in combination with an inter- 
active visualization of the restorations computed at each 
iteration. Unfortunately, visualization of CMB maps does 
not provide an accurate assessment of the accuracy of the 
computed solution. 

In the case of Tikhonov regularization, choosing a 
regularization parameter A is also a non-trivial issue. 
However, in comparison with stopping criteria for itera- 
tive methods, much more work has been done in this area 
llEngl et al. Il2n0nt iHanseTlliggTl: I\^^l20n2^ . The gen- 
eralized cross validation (GCV) method is probably the 
most well-known scheme for choosing a value for A. In 
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FWHM (arc.min.) 


Wiener 


Tikhonov 




rrms (%) 


rrms (%) 


a 1 a 


A 


10 


31.29 ±0.05 


30.88 ±0.06 


1.002 ±0.002 


0.73 ±0.01 


14 


33.62 ± 0.06 


33.07 ±0.06 


1.001 ±0.002 


0.75 ± 0.01 


23 


38.38 ± 0.08 


37.45 ± 0.09 


1.001 ±0.002 


0.79 ± 0.02 


33 


43.29 ± 0.09 


41.73 ±0.10 


1.001 ±0.002 


0.79 ± 0.03 



Table 1. Results of Tikhonov (with reflexive boundary conditions and discrete Laplacian for L) and Wiener deblurring 
of a Gaussian random field whose statistical properties are similar to those expected of the CMB sky observed with 
four channels of PLANCK-\^\ for beams with circular symmetry. The field has been contaminated with 100 different 
realizations of a white noise process [SjN = 2). The dimensions of the map are 340 x 340 pixels, which corresponds to 
a sky area of about 20° x 20°. The relative root mean square (rrms) is defined as the ratio of the residual root mean 
square (rms) to the rms of the true signal. The third column shows the ratio of the noise standard deviation estimate 
a defined by HA.4I) to the true noise standard deviation a. For the Tikhonov method, mean values and dispersions of 
the GCV estimates of A are also shown. 



FWHM (arc.min.) 


Wiener 


Tikhonov 




rrms (%) 


rrms (%) 


a/a 


A 


10 


30.97 ±0.72 


30.66 ±0.66 


0.96 ±0.02 


0.74 ± 0.04 


14 


33.24 ± 0.82 


32.80 ± 0.75 


0.96 ±0.02 


0.77 ±0.04 


23 


37.86 ± 1.00 


37.12 ±0.92 


0.95 ±0.03 


0.80 ±0.05 


33 


42.54 ± 1.18 


41.30 ± 1.07 


0.93 ±0.03 


0.81 ±0.05 



Table 2. As in Table with the only difference that a new Gaussian random field is generated for each simulation. 



Wiener Tikhonov 



FWHM (are.min.) 


rrms (%) 


rrms (%) 


a/a 


A 


10 


30.25 ±0.72 


29.96 ±0.68 


0.91 ±0.03 


0.73 ±0.03 


14 


32.29 ±0.81 


31.90 ±0.75 


0.91 ±0.03 


0.76 ±0.04 


23 


36.49 ± 0.98 


35.85 ±0.91 


0.90 ±0.04 


0.79 ±0.04 


33 


40.76 ± 1.15 


39.71 ± 1.04 


0.89 ± 0.04 


0.81 ±0.05 



Table 3. As in Table |21 but with a PSF of elliptical symmetry. The first column provides the FWHM along the major 
axis that is 1.3 times the FWHM along the minor axis. Axes of the PSF are parallel to the edges of the maps. 



this scheme, A is chosen to minimize the GCV function 



GCV(A) = 



Win 



[trace(J--H(A))/n]2' 



(21) 



where H. is the matrix that defines the estimator of if/, 
i.e., ytg ~ Hfxj n is the number of pixels in the 
image. For Tikhonov regularization 



n{x) 



(22) 



The development of the GCV method is based on the as- 
sumption that a good value of the regularization param- 
eter sh£uld_gredict missing data values (see Appendix 1X1 
and iGolub et al. '. .1979) . 

An additional argument for using Tikhonov regular- 
ization with GCV is that both have been well studied, 
and various authors have found this combination of regu- 



^Thompson et al. 


119911: iHanke & Hansen Ill993l: iHansen 1 


Il997: 


.Votrel..200^ 


). In Appendix lEl we describe how to 



efficiently compute both the minimum of the GCV func- 
tion (thus computing the regularization parameter) and 
the solution of the least squares problem (^0)1 . 

5. Numerical experiments 

We have used Monte Carlo simulations to check the re- 
liability and performance of the methodology presented 
in the previous sections. The simulations have been con- 
ducted under two different scenarios. We first fix the sky 
and generate different realizations of the noise process. 
Then, to account for the variability of the random field, 
we simulate different realizations of the random field and 
of the noise process. 

To keep the conditions of the experiment under con- 
trol but still relate it to the CMB problem considered in 
the next section, we have used Gaussian random fields 
characterized by a correlation function that has been ob- 
tained from the correlation function of the CMB via an 
approximation with an exponential function (for the de- 
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rmsw /rmsT 






FWHM (arc.min.) 


f3 = 1.0 


P = 0.5 


P = 0.25 


P = 0.15 


P = 0.10 


10 


1.01 


1.03 


1.10 


1.14 


1.16 


14 


1.02 


1.03 


1.10 


1.14 


1.16 


23 


1.02 


1.05 


1.11 


1.15 


1.17 


33 


1.04 


1.07 


1.16 


1.21 


1.23 



Table 4. Ratio rmsvy/rmsT for different values of the parameter /3 that increases the correlation length (see text). 



tails about the parameters used for CMB see Sect.|51). A 
dominating CMB component is a good approximation to 
real sky maps at least at medium and high Galactic lat- 
itudes where the effects of diffuse foregrounds from our 
Galaxy can be neglected Csee lMaino et al. II2OO2I and ref- 
erences therein). The PSFs used in the experiments are 
Gaussians with FWHM similar to that expected for the 
PLANCK-LFl instrument. However, since the exact form 
of the observing beam has not yet been determined with 
sufficient accuracy, we have tried two different scenarios: 
PSFs having circular and elliptical symmetry. For the lat- 
ter, the FWHM along the major axis is set to 1.3 times 
the FWHM along the minor axis. The axes of the beams 
are parallel to the edges of the maps. 

Since the random field is Gaussian and stationary and 
the noise is assumed white, classical Wiener filtering is 
expected to provide the smallest mean square error among 
linear filters. However, this filter requires knowledge of the 
spectrum of the unknown signal which is not available in 
practice. We use Wiener deblurring as a sort of benchmark 
to assess the performance of Tikhonov methods. 

Simulations were done using reflexive, periodic and 
zero boundary conditions and with two different types of 
penalty matrix L, the identity matrix and the discrete ap- 
proximation of the second derivative operator (note that 
realizations of a random field are smooth under m ild reg- 
ularity conditions on the correlation function, e.g., lAdler I 
However, since reflexive boundary conditions and 
the Laplacian operator have systematically provided bet- 
ter performance, only results concerning this combination 
are presented. In particular, we have found that the choice 
of the boundary conditions is a critical factor. In fact, not 
only have the results been systematically worse with zero 
and periodic boundary conditions, but with the latter we 
have also encountered stability problems in the estima- 
tion of the regularization parameter. These effects are the 
result of discontinuities introduced by zero and periodic 
boundary conditions; a problem that is much less imp or- 
tant for reffexive boundary conditions ijNg et aLlll999j) . 

The results of the simulations are shown in Tabs. ^ 
and 12 for circularly symmetric PSFs, and in Table 13 for 
PSFs with elliptical symmetry. Note that GCV estimates 
of the "optimal" value of A are stable with respect to both 
noise and different realizations of the random field. This is 
an important indication of the reliability of the methodol- 
ogy. It is well known, however, that GCV estimates may 
occasionally give very small values of A resulting in an 



under-smoothed solution, this problem can be easily cor- 
rected using a procedure defined in Appendix 1X1 

The same tables also show the mean value of the rel- 
ative root mean square (rrsm) error, defined as ||/ — 
/aI|2/||/||2- As expected, the error increases with the 
FWHM of the beam. The rrms error indicates that Wiener 
and Tikhonov deblurring (with refiexive boundary condi- 
tions and and discrete Laplacian) give comparable results. 
The advantage of the latter is that it does not require the 
spectrum of the unknown signal. That is, the smoothness 
constraint and the spectrum information provide similar 
results. (Note that since the Wiener filter has been imple- 
mented using the Fourier transform, only periodic bound- 
ary conditions were used for Wiener deblurring.) 

The third column in the tables compares estimates of 
the noise standard deviation obtained using ljA.4(l . a for- 
mula which arises naturally in the Tikhonov framework, 
to its true value a. These estimates are necessary to con- 
struct confidence intervals or test hypotheses. An over- 
estimate of cr indicates under-smoothing of the signal es- 
timate while an underestimate indicates over-smoothing. 
We see that Tikhonov provides reasonably good estimates 
of fJ. 

An important difference between Wiener and 
Tikhonov filters is that the former provides deblurred 
estimates assuming that / is a realization of a process 
with a particular spectrum, while the latter relies on 
the particular fixed realization of / on which the data 
are based. In particular, the Wiener fllter minimizes an 
error over all possible realizations of the signal while the 
Tikhonov filter assumes the signal is fixed. This means 
that while Wiener filtering relies on knowledge of the 
spectrum of the process to approximate the optimal 
filtering, Tikhonov uses the data to obtain an estimate of 
the signal's Fourier transform. Wiener filtering may thus 
give misleading results when the available realization 
of the signal is in the tail of the process distribution 
or when the particular realization of / happens to be 
somewhat unusual for the process, as it may happen 
when the assumed spectrum is incorrect. To illustrate 
this point, we repeated the simulations conducted for 
Table^but chose the initial sky realization from a process 
with successively larger correlation lengths; that is we 
incorrectly specified the spectrum required by Wiener 
filtering by multiplying the exponent of the correlation 
function by /?. Table 01 shows the ratio rms^y/rmsr of 
the rms error of the Wiener and Tikhonov deblurred 
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Fig. 1. Grayscale image of the simulated sky maps at 
30 GHz (see text). Each map contains 340 x 340 square 
pixels with side of 3.5' for a total area of about 20° x 20°. 

estimates for different values of f3. The results show that 
Tikhonov deblurring may provide better estimates when 
the spectrum is incorrectly specified. 

6. A CMB application 

We now present tests of the deblurring technique on simu- 
lated observations of the PLANCK-LF l. The region anal- 
ysed i s almost identical to the one in iBaccigahipi et al. I 
1(200 0^: it is a squared patch (340 x 340 pixels) with side 
of about 20°, centered at Z = 90°, & = 45° (Galactic co- 
ordinates). The latitude is high enough that CMB emis- 
sion dominates over foregro unds, assumed to be repre- 
sented by synchrotr on ((Haslam et al. I Il982l) and dust 
I Schlegel et al. iflQQSl) emission. We neglect contributions 
of point sources. The CMB model, in agreement with 
current experimental results l | de Bernardis et al. 1 1200 2: 
iHalverson et al. 1 120021: iLee et alT 200l|) . corresponds to 
a flat Friedmann-Robertson- Walker (FRW) metric with 
a cosmological constant (70% of the critical density), 
Hubble parameter today Hq ~ lOO/i km/sec/Mpc with 
h = 0.7 baryons at 5% and Cold Dark Matter (25% CDM), 
with a scale-invariant Gaussian initial spectrum of adia- 
batic density perturbations. 

The PLANCK-LFl instrument works at frequencies 30, 
44, 70, and 100 GHz. We assume nominal noise and angu- 



Fig. 2. Same as Fig.[T]but at 44 GHz (see text). 

lar resolution The maps are blurred through Gaussian 
PSF's with circular symmetry and appropriate FWHM's 
(i.e., w 33' at 30 GHz, w 23' at 44 GHz, « 14' at 70 GHz, 
w 10' at 100 GHz) and summed up together with simu- 
lated white noise with rms level as expected for the con- 
sidered channels. Since we choose to work with a pixel 
size of about 3.5 arcminutes, the noise rms are .042, .049, 
.042 and .043 mK in antenna temperature at 30, 44, 70, 
100 GHz, respectively. 

The original maps are shown in Figs. I1I4I together 
with the blurred, blurred plus noise, and deblurred ver- 
sions (Tikhonov method with reflexive boundary condi- 
tions and discrete Laplacian for L). The deblurred maps 
look reasonably good despite the high noise level. But, 
as expected, there is a clear loss of high frequencies, espe- 
cially for the lowest frequency maps. This loss, which is in- 
trinsic to any deblurring operation, is important given the 
high noise level. For example. Figs. EE show, respectively, 
the two-dimensional auto-correlation function of the origi- 
nal maps and their two-dimensional cross-correlation func- 
tions with the blurred (but noise-free) and the deblurred 
maps. These figures show that the deblurring operation 
provides good results for the lowest contour levels (those 
mainly determined by the lowest frequencies) and worse 
results for the highest contour levels ( those determined 
mainly by the highest frequencies). The same effect can 

* http:/ /astro. estec.esa.nl/SA-general/Projects/Planck 
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Fig. 3. Same as Fig.[T]but at 70 GHz (see text). 

also be seen in Fig. |H1 which shows a one-dimensional 
cross-section of the power-spectrum of the original 30 GHz 
and of the deblurred maps. 

The performance of the deblurring procedure can be 
also checked through the angular power spectrum, which, 
as usual, is defined by the expansion coefficients Cg of the 
two point correlation function in Legendre polynomials. 
Here t is the multipole associated to an angular scale of 
about 180/£ degrees. Consequently, our analysis applies 
to ^ ranging from t ~ 200, corresponding to the degree 
scale (larger angles are poorly probed due to the finite 
extension of our patches), up to multipoles corresponding 
to the instrumental resolution at each frequency. 

Estimates of the Cg coefficients for the original, 
blurred, blurred noisy, and deblurred maps are shown in 
Figs. 1911 21 Four CMB acoustic peaks are clearly seen in 
the original spectrum at ^ ~ 200,500,800, 1100. Since we 
are analysing a limited part of the sky, sample variance 
causes oscillations in the coefficients, with increasing am- 
plitude as the angular scale approaches the size of the 
patch, corresponding to the low multipole tail. 

As shown in the top panels of Figs. l9ll'21 the Cg are in- 
creasingly affected as the frequency decreases; the increas- 
ing tail at high multipoles is the effect of the instrument 
noise. In all four cases the deblurring procedure has two 
main effects. First, it restores amplitude and shape of the 
part of the spectrum not dominated by effects of instru- 



Fig. 4. Same as Fig.[T]but at 100 GHz (see text). 
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Fig. 5. Contour plot of the two-dimensional auto- 
correlation function of the original simulated CMB maps. 
Lag in pixels. 



ment noise and PSF. Second, as expected from Figs. l5l7l it 

reconstructs part of the power where the PSF causes a ma- 
jor degradation. In the 30 GHz case, the reconstruction is 
better in the multipole range 300 < £ < 400. Similarly, at 
44 GHz, anisotropy power is reconstructed up to £ ~ 500, 
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Fig. 6. Contour plot of the two-dimensional cross- 
correlation function of the original with the deblurred 
CMB maps. Lag in pixels. 
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Fig. 7. Contour plot of the two-dimensional cross- 
correlation function of the original with the blurred noise- 
free CMB maps. Lag in pixels. 
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Fig. 8. Onc-dimcnsional cross-section of the two- 
dimensional Power-spectra in two different frequency 
ranges of the original 30 GHz and the debiurred CMB 
maps (see text). Frequency is in Nyquist units. 
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Fig. 9. Angular power spectrum in different steps of the 
analysis at 30 GHz. In the top left panel the plot also 
shows the shape of the instrumental PSF in harmonic 
space (dotted line). 



and up to ^ ~ 700 at 70 GHz. Finally, at 100 GHz part of 
the original power is recovered up to £ ~ 800. 

Repeting the same simulations with PSFs of elliptical 
symmetry gives similar results. 

In Summary, in all four channels the deblurring pro- 
cedure was effective in recovering spectral properties of 
the maps. Although these results have been obtained un- 
der the restrictive assumptions of white noise, the char- 
acteristics of the deblurring procedure shown above are 
encouraging and deserve more attention and development 
in future work. 

A last comment regarding component separation. In 
principle, the separation can be obtained if maps of 
the same regions are available at different frequencies 



l)Maino et al. Il2002t iHobson et al. Ill998|) . In practice, this 
problem is technically difficult and is far from being 
solved. Its treatment is beyond the scope of this paper, 
but the important point is that maps must have the same 
spatial resolution in order to carry out the separation. 
This does not imply, however, that maps have to be de- 
blurred to a common resolution. One can also achieve a 
common resolution by blurring (smoothing). That is, we 
take Eq. lO and apply a frequency dependent smoothing 
operator S to both sides 

g = Sg = SS + Sz = j+~z. (23) 

Component separation can be done using the different es- 
timates of / obtained by smoothing the maps of different 
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Fig. 10. Same as Fig.Elbut at 44 GHz. 
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Fig. 11. Same as Fig. but at 70 GHz. 

frequencies. But without prior information on / the best 
Hnear estimates are simply the maps g. To include some 
prior information we can use smoothness constraints via 
Tikhonov regularization. The problem thus reduces to the 
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case we have considered before with the difference that 
this time the noise z is correlated due to the effect of 
S. Depending on the noise structure, one may be able 
to model correlated noise in Tikhonov regularization us- 
ing what is kno wn as mixed effect s mod els, some example s 
can be found in lMa fc Gul l|2002() ^ and lRobinsorTI l)l99l|) . 
If the noise is stationary, one can deconvolve in the fre- 
quency domain but again, when the signal spectrum is 
unknown we may want to regularize the problem using 
Tikhonov in the frequency domain. 



7. An example of deblurring effects on 
non-Gaussianity 

An important point in CMB research is the detection of 
non-Gaussianity. It is therefore of interest to check the 
effects of the deblurring procedure on the non-Gaussian 
characteristics of a map. 

Experimental results have not yet found any evidence 
of non-Gaussianity in CMB maps. It is even difficult to 
conduct realistic simulations of non-Gaussian CMB maps 
since it is not clear what type of non-Gaussian behavior 
one should look for. As an example, we consider the effects 
of deblurring on the marginal distribution of a particular 
homogeneous non-Gaussian random field whose marginal 
distribution is slightly different from a Gaussian (skew- 
ness ~ 0.22 and kurtosis ~ -0.11 - see Fig. The 
field has the same autocorrelation function, PSF (with 
circular symmetry) and S/N of the Gaussian fields simu- 
lated in Sect. 13 Realizat ions of this field ar e simulated via 
the method presented in I Vio et al. I l|200l|) . Skewness and 
kurtosis coefficients (normalized to be zero for a Gaussian 
distribution) are calculated for each simulation. 

Table El presents the average change, from the orig- 
inal values, in skewness and kurtosis over 100 simula- 
tions. As expected, noise and blurring have the effect 
of "Gaussianizing" the random fields (i.e., the estimated 
skewness and kurtosis are closer to zero than the origi- 
nal values), with obvious consequences on the response of 
any non-Gaussianity test. Deblurring noticeably improves 
estimates. 

These results, together with the power spectrum recon- 
struction discussed before, provide a further indication of 
the usefulness of deblurring in the analysis of CMB maps. 

8. An example of deblurring with spatially varying 
PSF 

Although the PSF of PLANCICs instrument is designed 
to be spatially invariant, this condition may change after 
launch, and there may be other experiments that require 
more general deblurring methods that can be used with 
spatially varying PSF. 

In principle, a real space approach to restoring an im- 
age degraded by a spatially variant PSF can be devel- 
oped, but in practice the methods are quite difficult to 



Fig. 12. Same as Fig.Elbut at 100 GHz. 



http://www.stat.purdue.edu/~chong 
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map 


A skewness 


A kurtosis 


30 GHz, blurred 

30 GHz, blurred noisy 

30 GHz, deblurred 


-0.052 ±0.058 
-0.098 ± 0.029 
-0.035 ± 0.024 


±0.022 ±0.123 
±0.055 ± 0.054 
±0.032 ± 0.054 


44 GHz, blurred 

44 GHz, blurred noisy 

44 GHz, deblurred 


-0.038 ± 0.040 
-0.088 ±0.017 
-0.030 ±0.019 


±0.021 ±0.087 
±0.054 ± 0.031 
±0.025 ± 0.042 


70 GHz, blurred 

70 GHz, blurred noisy 

70 GHz, deblurred 


-0.024 ± 0.024 
-0.078 ± 0.008 
-0.024 ±0.014 


±0.016 ±0.053 
±0.051 ±0.018 
±0.020 ± 0.032 


100 GHz, blurred 

100 GHz, blurred noisy 

100 GHz, deblurred 


-0.018 ±0.017 
-0.074 ± 0.008 
-0.022 ± 0.013 


±0.012 ±0.038 
±0.049 ± 0.020 
±0.018 ±0.027 



Table 5. Mean and standard deviation of the difference between the skewness and kurtosis coefficients of the blurred, 
blurrcd±noise, and deblurred non-Gaussian maps with respect to corresponding values (respectively w 0.22 and 
~ —0.11) of the original map (see text). The results are based on 100 simulations with S/N — 2. 



implement. Indeed, efficient impl ementations are currentl y 
available for iterative methods i Nagv & 0'Learv1ll997l) . 
but direct algorithms similar to those presented in previ- 
ous sections for invariant PSFs have not been developed 
yet. We are presently working on this problem. To illus- 
trate the importance of pursuing this work, we present 
some simulation results that use an iterative method to 
restore an image blurred by a spatially variant PSF. 

As claimed in Sect. 14. ll one of the most serious difhcul- 
ties with iterative methods is in the choice of an appropri- 
ate stopping criteria. Though methods such as the discrep- 
ancy principle, L-curve and GCV can be used, our success 
with these approaches on CMB maps has been marginal. 
Without an efficient method for choosing a regularization 
parameter, or for determining an appropriate stopping cri- 
teria, iterative methods may not be the ideal choice for 
CMB maps. However, because matrix-vector multiplica- 
tions with H and (the most expensive part in the 
iteration proc edure) can be done effi cientlv for spatially 
varying PSF llNagv fc 0'Le a,rvl l19C)7^ . iterative methods 
do give us a means of determining if better restorations 
can be obtained with a spatially varying model. 

To construct a test example, we have used a PSF 
that changes linearly across the map as shown in Fig. 
Deblurring was do ne using a stan dard conjugate gradient 
iterative method iHansenJllMl), 

stopping the iteration 
when the computed restoration is closest, in a mean square 
sense, to the true image. The matrix H modelling the spa- 
tially varying PSF is constructed by assuming that it is 
approximately spatially invariant in small regions. If the 
corresponding, spatially invariant matrix defined by the 
PSF in the ith region is Hi, then the spatially varying 
matrix H is defined as 



H 



D,H, 



where the matrices Di are nonnegative diagonal and 
= I, where I is the identity matrix. Note that using 



this H is equivalent to interpolation to match the indi- 
vidual PSFs across the different patches. For example, for 
piecewise constant interpolation the jth diagonal entry of 
Di is 1 if the jth pixel is in region i, and otherwise. 
Efhcient matrix vector multiplications with H and 
exploit the sparsity of £),j and the spatially invariant struc- 
ture of Hi. The implementation details are tedious to de- 
scribe; we only mention here that the basic idea is related 
to overlap-add and overlap-sa ve convolution methods, and 
refer the interested reader to llNagy fc 0'Lear v"l997^ for 
a description of the algorithms, and to l|Lee et al. Il2002|) 
for a Matlab implementation. 

The matrix H in our simulations was constructed 
through piecewise constant interpolation of 121 PSFs uni- 
formly distributed across the map. The computed spa- 
tially varying restoration and the corresponding spatially 
invariant restoration are shown in Fig. 1151 We have used 
a signal to noise ratio of S/N = 4 because the iterative 
procedure is more sensitive to noise than direct methods. 
Despite these limitations. Fig. 1141 shows the advantage of 
deblurring with a spatially varying PSF. 



9. Discussion and conclusions 

We have considered Tikhonov regularization for deblur- 
ring CMB maps in real space. Although more demanding 
from the computational point of view, this approach per- 
mits the development of algorithms that are more flexible 
and robust than those based on frequency-space methods. 
Furthermore, as shown in Fig. 1161 the computational cost 
can be significantly reduced by carefully implementing the 
algorithms to take advantage of the characteristic struc- 
ture of the matrices involved. 

We have applied the Tikhonov methodology to simu- 
lated skies at typical CMB frequencies. We considered test 
signals with known statistics, as well as realistic simula- 
tions of the CMB sky contaminated by noise whose rms is 
that expected for the low frequency instrument aboard the 
PLANCK satellite. This case is particularly interesting for 
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Fig. 13. Probability density function used in the simula- 
tion of the non-Gaussian random fields (see text). It cor- 
responds to a SB distr ibution belongin g to the Johnson's 
parametric family fsee lVio et al. l[r994l) . 



application of a deblurring procedure, as the instrument 
observes the sky at 30, 44, 70 and 100 GHz with very 
different PSFs of resolution 33, 22, 14, 10 arcminutes. 

We analysed the effects of the deblurring procedure by 
studying different characteristics of the restored image. 
Contour plots of the two-dimensional cross-correlation 
functions of the original and the deblurred maps show 
that the algorithm effectively improves the resolution, es- 
pecially that of the worst resolution channels at 30 and 
44 GHz. The same effect can be seen in the multipole co- 
efficients of the angular power spectra; the original power 
on angular scales hidden by the instrument's PSF is re- 
covered on a significant range of multipoles. We also per- 
formed an example of skewness and kurtosis recovery by 
the deblurring procedure. Instrumental noise and PSFs 
generally have the effect of pushing skewness and kurto- 
sis toward their Gaussian values; deblurring brings these 
values closer to their true non-Gaussian ones. 

On the basis of these promising results, we plan to 
develop deblurring algorithms based on the methods we 
have presented. Since Satellite CMB experiments are able 
to perform all-sky observations, a robust deblurring pro- 
cedure which is able to work on the whole sphere is con- 
ceivable. We also plan to test map based component sep- 
aration techniques, requiring multifrequency data of the 
same resolution, on deblurred maps. 



Appendix A: Automatic choice of the 
regularization parameter 

We briefly describe some methods to estimate the smooth- 
ing parameter A, which is an essential ingredient in 
Tikhonov regulariza tion . For further discussions of th is 
topic see iHansen. ((1997|) : iTenoTiol (jio'OL') : .Vogeil l|2002() . 
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Fig. 14. Some of the PSF's used in the experiment of the 
variant PSF described in the text. The couple of numbers 
in the header of each panel provides the coordinates of 
the pixels to which the displayed PSF corresponds. Each 
PSF is given by a two-dimensional Gaussian function with 
FWHM = 33 arcmin along the major axis and FWHM = 
10 arcmin along the minor axis. The size of each panel is 
23 X 23 pixels. 



Since the data g are noisy observations of Hf, it seems 
reasonable to choose a value of A that minimizes the pre- 
dictive mean square error (PMSE), 



PMSE(A) 



1 



Hf - Hf, 



(A.l) 



But since / is unknown, we minimize instead a cross- 
validation (CV) estimate of ljA.l|) obtained by plugging 
in ljA.l|) the data as a proxy for Hf and a leave-one-out 
estimate for H f ^. The result can be written as 



CV(A) = 



n ^ xi-n 



gi,\ 

(A) 



(A.2) 



where Ti.ii{X) are the diagonal elements of the "hat" ma- 
trix H22() that maps g into — Ti-f^. The CV estimate 
of A is the value that minimizes CV(A). 

The generalized cross-validation (GCV) function is a 
"smoothed" version of CV in which the diagonal elements 
of ?<(A) are replaced by their average 



GCV(A) 



9 dx Wl/n 



(1 - trace 7i(A)/n)2' 



(A.3) 



Unlike cross-validation, GCV is invariant under orthogo- 
nal transformations of the data. 

A slight modification of GCV provides and estimate 

of the noise variance 



9x 111 



72 — trace ?t(A) 



(A.4) 
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Fig. 15. a) original map; b) map obtained by blurring 
the map in the previous panel with a PSF changing uni- 
formly its orientation across the frame. Some of the PSF's 
are shown in Fig. (|14|l . White noise has been added to 
the map {S/N = 4); c) map obtained by deblurring with 
the spatial variant method explained in the text. A set 
of 121 PSF's uniformly distributed across the map have 
been used; d) map obtained by deblurring with a spatially 
invariant PSF. The PSF used is that shown in uppermost- 
left panel in Fig. H14|l . The standard deviation of the dif- 
ference between the deblurred and the true maps is 0.41 
for the spatially varying PSF and 0.43 for the spatially 
invariant one. 




600 700 
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Fig. 16. CPU time (sec.) required by the full Tikhonov 
deblurring procedure of iV x TV maps as a function of N. 
The PSF is the one used for the simulations of CMB maps 
at 30 GHz (circularly symmetric Gaussian FWHM 9.5 
pixels). Similar results are obtained for elliptic PSF with 
FWHM = 9.5 and FWHM = 2.8 pixels along the major 
and minor axes, respectively, and whose orientations are 
parallel to the sides of the map. Experiments have been 
conducted with Matlab 6 on a Pentium IV - 1500 GHz 
processor in a Windows 2000 operating system. 



minimum, it may lead to very small values of A (under- 
smoothing) . A lower limit on A that controls undersmooth- 
ing can be achieved by multiplying the trace term in IIA.3 [l 
by a constant fc > 1 l|Friedman fc Silverman lll989t ICu I 
I2OO2I) . 

The L-curve has not been studied as much as GCV 
but some studies seem to indicate that L-curve estimates 
have the same good properties of GCV a nd, in addition , 
may be more robust to cor r elated noise ijHansen Ill997l) . 
Note, however, that IVogel I l)l996|) has pointed out some 
convergence problems with L-curve estimates. 



This is just a normalized residual sum of squares where 
the effective degrees of freedom is determined by the trace 

of?t(A). 

The L-curve method l|Hansen Ill997j) is another way to 
determine a balance between goodness of fit and rough- 
ness. As A > increases, the points 



C{\)^{\\g-Hf 



\Lf,\\) 



(A.5) 



define a convex curve on the plane that in the log-scale 
resembles the letter "L" . The selection of A corresponds to 
the "corner" value where C(A) has the highest curvature. 

In general, estimates of A based on cross-validation are 
robust to small deviations from the homogeneous vari- 
ance and Gaussian assumptions and, for large samples, 
converge to the optimal minimizcr of the PMSE ( Wahba I 
Il990.) . When the GCV function is almost flat around its 



Appendix B: Efficient implementation of the 
Tikhonov algorithms 

B.l. Exploiting the structure of H and L 

The efficiency of our approach to computing a minimum 
of the GCV function l|21|l , and to solving the least squares 
problem (|2()(l . is based on exploiting structure of the ma- 
trices H and L. We assume that i is a structured matrix, 
such as given in Eq. ((T^ . 

Fast algorithms for certain structured matrices arising 
in image deblurring are well known. For example, when 
using periodic boundary conditions with a spatially invari- 
ant blur, the matrices H and L are BCCB, and therefore 
have the spectral factorizations 



H = .f*s:f. 



(B.l) 
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Here, T ^ Fn(E)Fm with "®" the Kronecker product and 
the onc-dimensional Fourier matrix that is a complex, 
unitary, and symmetric matrix whose elements are given 

by 

N and M are, respectively, the number of rows and 
columns of the image, ^* is the complex conjugate trans- 
pose of and S and A are diagonal matrices contain- 
ing the eigenvalues of H and X, respectively. Note that 
\i L = I, then A = J. The eigenvalues can be obtained 
by computing a two-dimensional fast Fourier transform 
(FFT) of the first columns of H and L, at a cost of 
0{N^ log A^), assuming the blurred image contains N x N 
pixels. Recently, however, many other efficient methods 
have been proposed that are suited to deal with a large 
variety of situations. 



B.2. Symmetric PSF 

If the spatially invariant PSF h{x,y) is symmetric, but 
not necessarily separable, i.e., h(x,y) — h{—x,y) — 
h(x,—y) — h{—x,—y), it happens that also the matrix 
H is symmetric . In this situation, it is possible to show 
liNg et a1. IIT99I that, under reflexive boundary condi- 
tions, the matrices H and L have the spectral factoriza- 
tions 



H 



C^AC. 



(B.3) 



where C is the orthogonal two-dimensional discrete cosine 
transform (DCT) matrix, and S and A are diagonal ma- 
trices containing the eigenvalues of H and L, respectively. 
In this case, the eigenvalues of H are given by 



[CHei] 



(B.4) 



where ef 



[ 1 • • • ] . Note that Hei is the first col- 
umn of H, which can be constructed from the PSF, and 
that multiplication by C can be done in 0{N'^ log N) op- 
erations using fast DCT algorithms. Computing the eigen- 
values Si of L is done in a similar manner. 

To efficiently compute the regularization parameter, 
we first replace H and L with their spectral factorizations 
in Eq. H21|) . and simplify to obtain 



GCV(A) 



/ 1:- 



(B.5) 

where, in the case of reflexive boundary conditions, g — 
Cg, and n — N"^ is the number of pixels in the image. We 
can now use standard minimization algorithms, such as 
Newton's method, to find the value of A that minimizes 
the scalar value function (|B.5|I . In the computations re- 
ported in this paper, we used Matlab's f minbnd function, 
which is based on Golden Section search and parabolic 
interpolation. 



We can also use the spectral factorization to efficiently 
solve the least squares problem (1^ . Because C is an or- 
thogonal matrix, H2U|) is equivalent to 



AA 



(B.6) 



where / = Cf and g = Cg (that is, / and g are, respec- 
tively, DCTs of / and g). Because S and A are diagonal 
matrices, this least squares problem can be solved, using 
a sequence of Givens rotations, with only O(iV^) opera- 
tions (Hansen 1997). In fact, strategic Givens rotations 
permits to change the structure of matrix 



AA 



( indicates a non-zero element) to 



(B.7) 



(B.J 



V / 

which is a form more amenable for an efficient solution. 

To summarize, after A is computed by finding the mini- 
mum of the GCV function, the Tikhonov solution (|12|l can 
be computed as follows: 

• Construct the first column of H from the PSF 

• Use a fast DCT algorithm to compute S. 

• Use a fast DCT algorithm to compute A. 

• Solve the least squares problem (|B.6|I to compute /. 

• Use a fast inverse DCT algorithm to compute / from 
/. 

The total cost of this approach is 0{N^ \ogN), and stor- 
age requirements are only 0{N'^) (i.e., the storage re- 
quired for anNxN image). Specific implementations used 
for the exper iments are availa ble in the Matlab package 
RestoreTools l|Lee et al. 11200^ ^ 

We remark that the efficient implementation described 
above assumes a symmetric H and reflexive boundary 
conditions. Of course a similar approach can be imple- 
mented for periodic boundary conditions, using FFTs in 
place of DCTs. However, if we want to use reflexive bound- 
ary conditions, and H is not symmetric, then other meth- 
ods should be considered. 



^ Available 
RestoreTools 



at |http://www.mathcs.emory.edu/~nagy/| 
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B.3. Separable PSF 

If the PSF is separable, then H can be decomposed into 
a Kronecker product, 



H = A(^B 



/ aiiB ai2B 

0-2lB 022-8 



ainB ^ 



(B.9) 



^B J 



The special block structure of Kronecker pro ducts can 
be exp loited is several ways ^ain 1989..; Kamm ^ Nasrv I 
Il998a|^ . In particular, assuming the images are N x N ar- 
rays of pixel values, then the following properties hold: 

1. The first important property is that we need only store 
the N X N matrices A and B, and do not need to 
construct the N'^ x N'^ matrix H explicitly. 

2. The X hnear system 

{A<^B)f = g (B.IO) 

is equivalent to the N x N matrix equation 

AFB^ = G, (B.U) 

where the N"^ x 1 vectors / and g are obtained through 
a lexicographical ordering of the N x N arrays F and 
G. 

3. The singular value decomposition (SVD) is a tool used 
for an alyzing and solving ill-posed problems ( Hansen I 
I1997I) . It is usually too expensive for large scale prob- 
lems, such as image deblurring. However, for Kronecker 
product structures, the SVD can be used efficiently. To 
see this, suppose 



and B = Ub'SbVl 



(B.12) 



A = Ua-SaVi 

are the SVDs of A and B. Then 

A®B = (t7a®t7fc)(S„0Sfa)(yf (B.13) 

is the SVD oi H = A(E)B.ln particular, the SVD of a 
Kronecker product can be computed at a cost of only 
0{N^), rather than 0{N^) if working directly with H. 

Using these properties, it looks like the SVD can be used 
in place of the the spectral factorization of H. However, 
we run into difhculty if the regularization operator, L, is 
not the identity matrix. To see this, assume that the SVD 
of H is given by 



H = UH'SVi,, 



(B.U) 



where Uh ^Ua®Ub,'^^ and Vh = Va®Vb- 

In addition, assume L has the SVD 



(B.15) 



by the above properties of Kronecker products, we see 
that the cost of transforming (|21|l into l|B.5|l and H20I) into 
(|B.6|I is 0{N^). This is slightly more expensive than the 
0{N^ log N) cost when using fast transforms, but it is still 
a reasonably efficient approach. The storage requirements 
remain 0{N^). 

A problem arises, though, when trying to use a differ- 
entiation operator for L because, in general, Vh 7^ Vl- 
If these two orthogonal matrices are not equal, then H20|l 
cannot be efficiently transformed into ljB.6|l . Moreover, 
H21|l cannot be efficiently transformed into ljB.51) . and 
so evaluation of the GCV function, and hence computa- 
tion of its minimum, becomes much more expensive. In 
this difficult situation, a hybrid iterative/direct method 
((Kilmer fc 0'Lear'v1ll999l) may be appropriate. 

Finally, we remark that in general there is no computa- 
tional difference between separable spatially variant PSFs 
and separable spatially invariant blurs. That is, efficient 
direct methods can be used if L = J, but not when usin g 
a differentiation operator for L l(Kamm fc Nagv (ll998b|) . 
Moreover, in the difficult cases when the direct factoriza- 
tion approach cannot be used, iterative and hybrid meth- 
ods can be implemented very efficiently f or both spatially 



invariant and spatialh 


r varying blurs 


Il99fit 


iNa.gv fc O'T.earv 


Il997lll998l). 



where the matrices Uh , Vh , U l , and V l are orthogonal, 
and S and A are diagonal, li L = I, then A = /, and 
we can take Ul = Vl = Vh- Now Eqs. (|B.5|I and (|B.6p 
can be used with / = Vjjf and g — U^g. In this case, where: 



B.4. Non-separable PSF 

In this section we have so far described efficient implemen- 
tations of direct methods for the following situations: 

1. Spatially invariant PSF, with periodic boundary con- 
ditions. In this case FFTs are used. 

2. Spatially invariant and symmetric PSF, with reflex- 
ive boundary conditions. In this case fast cosine trans- 
forms are used. 

3. Separable PSF, invariant or variant. In this case ef- 
ficient use of the Kronecker product structure is ex- 
ploited. 

If the image deblurring problem does not fit into one of 
these categories, then some other approach should be used. 
One possible option is to simplify, or approximate, the 
problem with one that does fit into one of the above cate- 
gories. For example, a symmetric approximation, such as 
{H + H'^)/2, can be used with reflexive boundary con- 
ditions so that the fast cosine transform approach can be 
used. Other, recently developed approaches, use more so- 
phisticated approximation techniques. 

In particular the Kro necker approximation method 
llKamm fc Nagv I Il998albl) is a technique that deserves 
some attention. The main idea of this method is to ap- 
proximate the non-separable matrix H with a separable 
matrix. To see how this can be done, suppose P is an rt x n 
image of the PSF, with pij denoting the center (origin) or 
the PSF. If the PSF is separable, then we can write 



P = ab^ and H = A®B , 



(B.16) 
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Zero boundary conditions imply A and B are Toeplitz 
matrices defined by a and b, respectively. That is, 



/a, ■■■ ai 



and 



V 



B 



bn 



V 



ai 

a„ • ■ • ai J 
\ 

bi 

bj 



(B.17) 



(B.18) 



Reflexive boundary conditions imply A and B are 
Toeplitz-plus-Hankel matrices defined by a and b, re- 
spectively. That is, 



/ Gi ■ ■ ■ ai 



0.1 



\ an ■ ■ ■ ai J \ ai ■ ■ ■ a^-i 

(B.19) 



and 



B 



(b, ■■■ bi 



\ ( bj+i 



bi 



bn ■■■ bjj \ 



bi 



bi ■ ■■ bj-i 
(B.20) 



If the PSF is not separable, then we could compute a 
rank-one approximation of P, and construct A and B as 
described above. That is, 



P«ab^ ^ H^A(E)B. 



(B.21) 



With a proper weighting applied to P, one obtains an 
optimal approximation of the form 



min Hi? - A (g) BI 



(B.22) 



where || • is the Frobenius norm and the minimiza- 
tion is done over all Kronecker products A(E) B. This ap- 
proximation can be computed in 0(n^) operations; see 
iKamm fc Naev I l|l998 a b) for further details. 
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